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A semi-Lagrangian method for advection equation with hy- 
brid cubic-rational interpolation is introduced. In the present 
method, the spatial profile of physical quantities is interpo- 
lated with a combination of a cubic and a rational function. 
For achieving both high accuracy and convexity preserving 
of solution, the two functions are mixed in the optimal ra- 
tio which is given theoretically. Accuracy and validity of this 
method is demonstrated with some numerical experiments. 

Key Words: Numerical method, Advection, Semi- 
Lagrangian method, Interpolation, Cubic function, Rational 
function, Convexity preserving. 



I. INTRODUCTION 

Basically, the CIP method ||,g] is a numerical method 
for an advection equation, 

df , df 

Tt +u Tx =Q 

whose solution is expressed as 

f(x,t + At) = f(X(x,t),t), (1) 

where X is the trajectory of fluid particle, which is lo- 
cated at x at the time t + At, 



f(x, t + At) = f{x - u(x, t)At, t), 



(3) 



X{x,t) = x + 



t-At 



u(X(x, t), T)d/T. 



(2) 



In a semi-Lagrangian scheme like the CIP, the solution 
(|l|) is solved as an interpolation problem |||l[ . In the CIP 
scheme, an Hermite cubic expansion function is used to 
interpolate / at the time t. The quantity and its first 
spatial derivative defined at each grid points are updated 
as to obey, respectively, eq. (|I|) and its spatial derivative, 
i.e., 

df{x,t + At) dX{x,t) df(X(x,i),t) 



dx 



Ox 



dX 



Generally, the integration in eq. (0) is solved by assuming 
that u is locally constant as 

/t-At 
(It 

— x — u(x, t)At. 

With this, the solutions (pi) and (0) are expressed, re- 
spectively, as 



df(x, t + At) _ du{x,t) df_(x -u(x, t) At, t) 



Ox 



dx 



dX 



(4) 



In the last decade, various kinds of extension and 
improvement have been adopted to this method. In 
1990, Yabe et al extended this to multidimensions with- 
out time-splitting technique by employing a multidimen- 
sional cubic expansion function M. In 1991, Kondoh ex- 
tended the ID CIP to a 5th-order advection method and, 
furthermore, proposed a solver for parabolic equations by 
extending the basic concept of the CIP || . In the same 
year, Aoki and Yabe improved the multidimensional one 
by modifying the multidimensional expansion function 
and proposed two alternative formulae ||[i]]. In 1994, 
Kondoh extended his approach to a multidimensional 
parabolic equation and general hyperbolic equations jB|. 
In 1995, Ida and Yabe proposed an implicit version of 
the CIP ||. This method is CFL free and can be solved 
directly with a marching procedure although it is a 3rd- 
order method. In the same year, Utsumi extended the 
CIP to a solver for the Euler equations of fluid flow with- 
out finite-difference technique by employing differential- 
algebraic and Lagrangian-like concepts pOJ . In 1996, 
Xiao et al proposed a convexity preserving method for 
the advection equation by replacing the cubic function 
in the CIP scheme with a cubic-rational function PtL[12| . 
In the same year, Ida proposed a high accurate solver 
for free-surface flow problem by coupling the CIP with 
newly proposed extrapolation scheme p3|-[l5[ . With this 
method, the density discontinuity at material interface 
is solved without any numerical dissipation across the 
interface. In 1997, by extending the Kondoh's approach 
B , Aoki proposed high accurate solver for wave equation, 
full Euler equations and others |l(| . In 1999, Tanaka et al 
proposed an exactly conservative solver for the continuity 
equation in non-conservative form by additionally using 
the local mass of fluid as a dependent variable 0,0 ■ 

In this paper we discuss on the rational method pro- 
posed by Xiao et al. As proved theoretically in ref. 11, the 
rational interpolation method suppresses the numerical 
oscillation which tends to appear in high-order solution. 
However, this method sometime provides more diffusive 
result than that with the classical cubic interpolation. 
We try to improve its accuracy, without any loss of the 
convexity preserving property, by mixing with the cubic 
interpolation function in the optimal ratio. In Sec. 2, the 



conventional rational method is briefly reviewed and, in 
Sec. 3, the optimal mixing ratio is given theoretically. In 
Sec. 4, results of some numerical experiments are shown 
for demonstrating the accuracy and the validity of the 
present method. 



II. CONVEXITY PRESERVING 
SEMI-LAGRANGIAN METHOD 



( CR(0,l) = fi, 

CR{h,l) = fi+i, 

d x CR(0,l) = d h 

[d x CR(h,l) = d i+1 . 



The rational formula is adapted in a cell where the data 
is convex or concave, i.e., 

di > Si > di -Li 



In the rational method [0, the following cubic- 
rational interpolation function is used: 



CR&i) 



fi + AU + A2 % e + A3 % e 
1 + 7^ 



(5) 



where 



Ah = di + faB h 
A2i = SaBi + (Si - di)/h - A3ih, 
A3 t = [di - S t + (di+i - 50(1 + lB t h)]/h 2 
n i Si dj 

Di = 



Si 



di+i — 5 
(Ji+i - fi)/h 
—UiAt, 



-1 \/h 



fi and di are a physical quantity and its first spatial 
derivative, respectively, Ui is the particle velocity at Xi 
assumed as negative here, At is the time interval, h is 
the grid width assumed as uniform in this paper for sim- 
plicity and 7 is a parameter for switching the form of the 
interpolation function. In the case of 7 = 1, eq. rtq) is 
reduced to a rational formula of 



CRfol) 



fi + Rig + R2 % e 
1 + Bit 



(6) 



where 



Rl t — di + fiBi, 
R2i = SiBi + (Si 



di)/h. 



Unlike rational functions used in data interpolation tech- 
nique (See Ref. [|l9| for example), the above formula is 
constructed not only with the quantity but also with its 
derivative. In the case of 7 = 0, on the contrary, eq. (0) 
is reduced to a cubic formula of 



CRfo 0)=fi + CUt + C2 t e + CM 3 



(7) 



where 



cu 


= di, 


C2i 


= -(2di + d i+1 -3Si)/h 


C3 t 


= (di + d l+l - 2S l )/h 2 . 



Those rational and cubic formulae satisfy a continuity 
condition at Xi and X%+i expressed as 



di < Si < d l+1 . 

This rational interpolation function preserves the convex- 
ity of solution. For the other data, interpolation function 
is switched to the formula of eq. (7\) which corresponds 
to the conventional one of CIP Q and provides purely 
3rd-order solution. 

With those interpolation functions, in this paper, we 
propose a hybrid method of making use only of their 
superior characteristics by mixing them optimally. In 
ref. 12, Xiao et al proposed an additional switching tech- 
nique shown as 



7 



1, for di ■ di+i < 0, 
0, otherwise. 



(8) 



This means that the rational function is applied only in 
the cell which includes a turning point of the gradient. 
While this procedure modifies the dissipation property of 
the rational method, this breaks the preserving of con- 
vexity as shown in Sec. 4. The optimal mixing technique 
being proposed in this paper would achieve improvement 
of accuracy without any loss of the convexity preserving 
property. 

For the convenience of the following discussion, we re- 
arrange eqs. (ph and (fjj) as 



R(k) = fi + dM + 



P?k 2 



Qi + (Pi-Qi)k' 



(9) 



and 



C(k) =fi + d t hk + (2P, - Qi)k 2 + (Qi - P t )k 3 , (10) 
respectively, where 

Pi = (Si - dt)h, 

Qi = (d i+ i - S t )h, 

and 

k = £/h = -mAt/h 
is the local Courant number and 

k e [0, 1] 

because of the CFL condition. 



III. HYBRID CUBIC-RATIONAL METHOD 
WITH THE OPTIMAL MIXING 

A. The optimal mixing of the two interpolation 
functions under the convexity-preserving condition 

We start from a combination of the rational and the 
cubic functions shown as 



F(k) = aR(k) + (l-a)C{k), 



(11) 



where a is a weighting parameter whose range is limited 
as a G [0,1]. For large a, i.e., a ~ 1, less oscillatory 
solution may be expected because the rational function 
becomes dominant and, for small a, i.e., a ~ 0, high- 
order solution may be expected because the cubic one be- 
comes dominant. By determining a properly, convexity- 
preserving high-accurate method may be produced. We 
discuss below how to determine the weighting parameter. 
While numerous proofs on some characteristics of the 
rational interpolation method have been done in Rcf . pi| , 
the roof of all of them is a property of the rational func- 
tion. The property is expressed as that, for the convex 
data, the interpolation function is convex between the 
given interval and, for the concave data, it is concave. 
Namely, 



d 2 R{k) 
dx 2 



<0, forfce[l,0] 



is true for the convex data of 



Pi < and Qi < 0, 



and 



d 2 R{k) 
dx 2 



>0, fork e [1,0] 



is true for the concave data of 

P % > and Qi > 0. 

For achieving both convexity preserving and better res- 
olution, we make the weighting parameter the minimum 
value which satisfies the above condition. The follow- 
ing discussion will be limited to the case where the data 
is convex or concave since the rational interpolation is 
adopted only in the case. For the other data, the cubic 
interpolation is adopted. 

The second spatial derivative of eq. (fTD) is shown as 



d 2 F{k) 
dx 2 



a 



« 2 Q? 



2P 2 Q 



h 2 [(Pi - Qi)k + Qi} 3 



.(l-a) — [2,{Q i -P i )k + 2P i -Q i }. (12) 



By introducing Hi defined as 



eq. (y_2|) is rewritten as 
d 2 F(k) _ 2 



H, 



dx 2 h 2 ^ l \ [{l-H^k + H.Y 

3(H t - l)fc + 2 - H, 



+(l-a)- 



H, 



For preserving convexity of solution, as mentioned above, 
the following condition must be satisfied: 



or 



d 2 F{k) 
dx 2 



d 2 F(k) 
dx 2 



< 0, for P t < and Q t < 0, 



> 0, for P t > and Q t > 0. 



This condition is expressed by an inequality of 

Hi 



d 2 F(k) /(2_ ( 



dx 2 



' [(1 - Hi)k + H % Y 



3( Hl -l)k + 2-H t ^ 
+ (1 - ") jj, > 0, 



(13) 



because the sigh of d 2 F(k)/dx 2 and Qi should be the 
same. Furthermore, the inequality (|l|) can be rewritten 
as 

aH 2 + (1 - a)G 3 i[3(Hi - l)fe + 2 - Hi] > 0, (14) 

where 

Gi = (1 - Hi)k + Hi, 

because of 

Hi > (15) 

(Remind that the signs of Pi and Qi are the same) and 

Gi>mm(l,Hi) > 0, for fee [0,1]. (16) 

Here we introduce Ki(k) which is defined as 

Ki(k) = G?[3(Hi-l)k + 2-Hi\. 

With this function, eq. ( fl4| ) is rewritten as 

aH? + (1 - a)Ki{k) > 0. (17) 

From the inequality (flq) and 

3(iJ. t - l)k + 2 - Hi > min(2 - Hi, 2H L - 1), 
forfc e [0,1], 

it is known that the inequality (117n is always true for 

1/2 < Hi < 2 

because of Ki(k) > 0. In this case, we can use a = 0, i.e., 
the purely cubic interpolation function which provides 
high-order solution. 



For a while, the following discussion will be limited to 
the case of 

Hi >2. 

The first and the second spatial derivatives of Ki(k) re- 
spect to k are 



dKjjk) 
dk 



= 6(l-H l ) 2 G 2 (l-k) 



(18) 



and 



d 2 K % (k) 
dk 2 



= 6{l-Hi) 2 Gi{2-3Gi), 



(19) 



respectively. From eq. (fill) with an inequality of 

d > min(l, Hi) = 1, for Hi > 2 and k G [0, 1], 
we know that 

d 2 K t (k) 



dk 2 



<0, for fee [0,1], 



which means that dKi(k) /dk is monotonicaly decreasing 
for k £ [0, 1]. Furthermore, from eq. (|l8|), we know that 



dKj(0) 

dk 



= 6(1 - HifHf > 



and 



dKjjl) 

dk 



0. 



Namely, dKi(k)/dk has non-negative value for k £ [0, 1]. 
Those results show that Ki(k) is monotonicaly increas- 
ing for k £ [0,1], and the minimum value of it is 

min(K t (k)) = Ki(0) = H?{2 - H t ) < 0. 

Thus, the inequality ( fi"7| ) is reduced to 

aHf + (1 - a)H?(2 - H t ) > 0. 

From this, we get 

HiiHi-2) 



a > 



Hi(Hi-2) + l 



For increasing the dominance of the cubic function as 
much as possible, we use 



= HjjHi - 2) 
Hi(Hi-2) + r 

The proof for the case of 

< H t < 1/2 

can be done easily as follows. With replacements of 



(20) 



Hi - 1/Li 
k->l-l, 



(21) 



the inequality (|lj) is rewritten as 

aL 2 + (1 - o)[(l - U)l + i,] 3 [3(L 4 -1)1 + 2- L t ] > 0. 

This inequality corresponds to that of (|14| ) only with re- 
placements of Li — » Hi and I — » k. Thus, the solution of 
this inequality is 

a > L ; T {Ll ;. 2 | - , forL«>2andie[0,l]. (22) 

From eq. (|2l|), it is known that the result (E2n is for 
< Hi < 1/2 and fc G [0, 1], From the same reason as 
mentioned in leading eq. (p0|), we use 



Li(Lj-2) 
Li(Li-2) + r 



(23) 



The above results on the determination of a is sum- 
marized as 



where 



MjjMj-2) 
M i (M i -2) + l' 



Mi = max[2, max(_ff i , L,-)] 
= max[2, max( — , — — )1. 



(24) 



(25) 



In the case of 1/2 < iJ, < 2 and i/^ < 0, M, and a are 
determined with eqs. (|25| ) and ( p4[) as 

Mi = 2, 
a = 0, 

and the interpolation function is reduced to pure cubic. 
For H t > 2, eq. @ is reduced to eq. @ and, for L % > 2, 
it is reduced to eq. (E3T). 



B. Summary of the formula 

With a little rearrangement, the formula derived in the 
last subsection is summarized as 

F(k) = f l + dM+{Gh + G2 l )k 2 , 
— — = d i + [Gl l — \-2G2i + (l-a){Qi-Di)\~, 



where 



GU = aPf/D,, 
G2i=(l-a)(2Pi-Di), 
D % = Q % + (Pi - Q t )k, 

M(M-2) 
M(M-2) + l' 

Mi = max[2, max( — , — ) 



with 

Pi = (Si - ck)h, 

Qi = {d-i+i — Si)h, 
Si = (fi+i - fi)/h. 

For the case of Ui > 0, we need replacements of 

i 4- 1 — ► i — 1, /i — * —h. 



where the superscript n shows the number of time step. 
In Fig. 1, we plot the error as a function of the grid width 
ft at n = 4000 with CFL = 0.2. From this result, it is 
known that the accuracy of the methods except for the 
conventional rational method are very similar each other 
on this problem, while that of the present one is inferior 
a few among them. 
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FIG. 1. Numerical error as a function of the grid width h in 
the results with the present, the CIP, the conventional ratio- 
nal and the modified rational (with the additional switching 
technique) methods. 



IV. NUMERICAL EXPERIMENTS 

For demonstrating the validity of the previous discus- 
sions and the accuracy of the present scheme, we show 
some numerical experiments in this section. 

The first example is the linear propagation of a sinu- 
soidal wave. The initial condition is 

f(i,0) = 0.5cos(4.7r/ii), 

where the grid width h is set in this example as 

h = 1/JV 

and N is the number of grid points. The velocity is set as 
u = 1 and is assumed as constant in space and time. The 
boundary condition at i = and N is periodic. With 
this example, we compare the accuracy of the present 
method with that of the three existing methods, i.e., the 
CIP, the conventional rational method and that with the 
additional switching (||). The last one is called below 
the modified rational method. Here we define numerical 



Error = 



N 

J2 |/™ — analytical value | 
»=i 

N 
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FIG. 2. Linear propagation of a triangular and a square 
waves. The figure shows results at n = 1,000 (the upper 
group) and n = 10, 000 (the lower group) with the present (1), 
the CIP (2), the conventional rational (3) and the modified 
rational (4) methods. 

Next we solve the linear propagation of a triangular 
and a square wave. The pulse widths of the waves are 
30/i and 26/i, respectively, and pulse heights are 1. In 
Fig. 2, we show results at n = 1, 000 (a) and 10,000 (b) 
with CFL = 0.2. In the results with the CIP method, 
overshoots and undershoots are obviously shown around 
the discontinuities of the square wave. The results with 
the conventional rational method are more diffusive than 
other ones. The results with the present and the modified 
rational method are very similar while the former one is 
more diffusive quite a little and the latter one has small 
over- and undershoots. 

With the same example of the square wave, we fur- 
ther discuss on the convexity preserving property of those 



methods. Here we introduce a variable p defined as 
f 1/2, if/f +1 -/f>0, 

p-Vi/2 = < o, if /p +1 - /f = o, 

[-1/2, if/f +1 -/f<0. 

By observing this variable, we can appreciate the con- 
vexity of solution. In Fig. 3, we show / and p at n = 150 
(a) and 10,000 (b). Analytically, p becomes 1/2 and 
— 1/2 only at the left and right discontinuities, respec- 
tively, and p = elsewhere, namely, only one positive 
and one negative regions should exist in the spatial pro- 
file of p if the convexity of / is preserved. In the results 
with the present and the conventional rational method, 
the two regions are clearly seen while the width of the 
regions is expanded by the numerical diffusion. In the 
results with the CIP and the modified rational methods, 
many numbers of leaps are shown in the profile of p. This 
means that the conventional CIP and the modified ratio- 
nal methods are not oscillation free. 
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FIG. 3. Results of the linear propagation of a square wave 
at n = 150 (the upper group) and n — 10, 000 (the lower 
group) with the present (1), the CIP (2), the conventional ra- 
tional (3) and the modified rational (4) methods. The circles 
and the solid lines show / and p, respectively. 





(b) Step =10,000 



FIG. 4. Results with the reduced mixing ratio (0.99, 0.95, 
0.9, 0.8 and 0.7 times values of the optimal one) at n — 150 
(the upper group) and n = 10,000 (the lower group). 

In order to demonstrate that the determination of a 
with eq. ( p4{ ) is the optimal, we show results of the same 
example with the present method and reduced a. In 
Fig. 4, we show results with 0.99, 0.95, 0.9, 0.8 and 0.7 
times smaller value of a than that is determined with 
eq. (p3) and CFL = 0.2. The unphysical oscillation is 
observed even in the results with 0.99 times value and 
the oscillation becomes stronger according as a decreases. 
The numerical oscillation at n = 10, 000 becomes weaker 
than that at n — 150 because of numerical diffusion, but 
still exists. The overshoot around the discontinuity of / 
can also be seen in those results and it becomes larger as 
a decreases. 

Finally we show results of an extreme example for 
demonstrating the robustness of the present method and 
a defect of the modified rational method. The initial 
profiles of / and u, the latter is assumed as constant in 
time, for this example are shown in Fig. 5. In the initial 
profile of /, there exist two steep gradients which are ap- 
propriately smooth but sufficiently sharp. In the velocity 
distribution, a similar gradient exists, and u = 1 in the 
left side of the gradient and u = 0.1 in the right side. 
Those mollified profiles are made as follows: First / and 
u are set as 



1, 



for 5 < i < 67, 
elsewhere 



and 
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FIG. 5. Initial condition of the extreme example. The up- 
per and the lower figures show / and u, respectively. 



In solving this example, we estimate the velocity gradi- 
ent appearing in eq. (gj) with a conventional second-order 
centered finite differencing as usually done in the CIP 
scheme JjJ. 
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FIG. 6. A typical result of the extreme example by using 
the present method. The figure shows results at n = 0, 140 
and 450 with CFL = 0.25. 



1. 
0.1, 



for i < 71, 

elsewhere. 



Next, those values are smoothed with a conventional 3- 
point smoother, 



(m+l) _ 



(m) 



(to) . (m) 



= (l-e) 9 r+e 
for m = 1,2,- --M, 

where g is / or u, m is the iteration number of the 
smoothing, M is the maximum of it and e is a pos- 
itive constant smaller than 1.0. For / and u, we use 
(M,e) = (2,0.05) and (2,0.1), respectively. The result- 
ing mollified values are used as initial values. Initial value 
of dj is set as 



0, 



(M) 



f (M) 



if/ (M) = 0or/ (M) = 1) 



fl_i)/2h, elsewhere. 



FIG. 7. Results when the right discontinuity just passes 
through the velocity gradient (the left group:n = 48 with 
CFL = 0.25. the right group: n = 24 with CFL = 0.5). In 
the results with the CIP and the modified rational methods, 
strong overshoot appears. 
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FIG. 8. The maximum (upper) and the minimum (lower) 
values in the numerical results of the extreme example as a 
function of the time step. The maximum and the minimum 
values calculated with the present and the conventional ratio- 
nal methods keep with 1 and 0, respectively, while those with 
the CIP and the modified rational method are oscillated. 

In the first stage of this problem, the mollified square 
wave of / propagates to right in u = 1. Then, after the 



wave covers the gradient of u, the waveform is compressed 
in horizontal direction and the pulse width is decreased. 
Finally, after passing through the velocity gradient, the 
width of the square becomes 1/10 of the initial one, i.e., 
~ 6ft, and the square propagates in u = 0.1 (See Fig. 6 
in which a typical result of this problem by the present 
method is shown). 

In Fig. 7, we show results at n = 48 with CFL = 0.25 
and at n — 24 with CFL = 0.5. While the results with 
the present and the conventional rational methods are 
smooth and oscillation free, strong overshoot is shown in 
those with the CIP and the modified rational methods. 
At the time where the gradient of / passes through the 
gradient of u, it is amplified strongly by the velocity gra- 
dient. Thus, if the cubic interpolation is adopted in this 
case, the overshoot should appear. In Fig. 8, we show 
the maximum and the minimum value of / as a func- 
tion of the number of time step. The maximums with 
the CIP and the modified rational methods have strong 
peak at n ~ 50, i.e., when the right gradient of / passes 
through the velocity gradient and, furthermore, the min- 
imum with the CIP has strong peak at n ~ 300, i.e., 
when the left gradient of / passes through there. This 
result shows that the additional switching (B) raises in- 
adequate adaptation of the cubic function. The present 
method does not have such a defect. In Fig. 9, the results 
at n = 550, i.e., those after the pulse has passed through 
the velocity gradient is shown. The result with the mod- 
ified rational method has weak overshoot and that with 
the conventional rational method is more diffusive than 
that with the present one. 

The above results prove that the present method has 
higher accuracy than that of the conventional rational 
method and is a convexity-preserving method. 
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FIG. 9. Results after the square have passed through the 
velocity gradient (n = 550 with CFL = 0.25). The right figure 
shows an enlarged one of a part of the left one. The result 
with the CIP is not shown in the right figure. Theoretically 
the height and the width of the projection should be 1 and 
about 6ft, respectively. 



V. CONCLUSION 

In this paper we proposed a hybrid semi-Lagrangian 
method with a cubic and a rational interpolation func- 
tions. The optimal ratio for mixing those functions was 
led theoretically. The present method has higher accu- 
racy than the conventional rational method and is os- 
cillation free. The numerical experiments curried out in 
the last section demonstrate the validity of the theoret- 
ical discussion in Sec. 3 and the accuracy of the present 
method. However, the results show a limit of the present 
method as well. Higher accuracy than that of the given 
results shown, for example, in Figs. 2 and 3 may no 
longer be expected when one use the cubic and the ratio- 
nal functions under the convexity preserving condition. 
Extension to a higher-order method may be needed for 
achieving higher resolution. 

The present method may be extended to multidimcn- 
sions and a conservative method by employing the ap- 
proaches discussed by Aoki |||7| and Tanaka et al |17|,[L8| , 
respectively. 



[1] T. Yabe and T. Aoki, A universal solver for hyper- 
bolic equations by cubic-polynomial interpolation I. One- 
dimensional solver, Comput. Phys. Commun. 60 (1991) 
pp.219-232. 

[2] T. Yabe, T. Ishikawa, P. Y. Wang, T. Aoki, Y. Kadota 
and F. Ikeda, A universal solver for hyperbolic equations 
by cubic-polynomial interpolation II. Two- and three- 
dimensional solvers, Comput. Phys. Commun. 66 (1991) 
pp. 233-242. 

[3] D. R. Durran, Numerical methods for wave equations in 
geophysical fluid dynamics, Springer, Sec. 6, p. 303. 

[4] T. Yabe, T. Ishikawa and Y. Kadota, A multidimensional 
cubic-interpolated pseudoparticle (CIP) method with- 
out time splitting technique for hyperbolic equations, F. 
Ikeda, J. Phys. Soc. Japan 59 (1990) pp.2301-2304. 

[5] Y. Kondoh, On thought analysis of numerical scheme for 
simulation using kernel optimum nearly-analytical dis- 
cretization (KOND) method, J. Phys. Soc. Japan. 60 
(1991) pp.2851-2861. 

[6] T. Aoki and T. Yabe, Multidimensional cubic interpola- 
tion for ICF hydrodynamics simulation, NIFS-82, 1991, 
pp. 1-9 (unpublished). 

[7] T. Aoki, Multi-dimensional advection of CIP (Cubic- 
Interpolated Propagation) scheme, CFD J. 4 (1995) 
pp. 279-292. 

[8] Y. Kondoh, Y. Hosaka and K. Ishii, Kernel optimum 
nearly-analytical discretization (KOND) algorithm ap- 
plied to parabolic and hyperbolic equations, Comput. 
Math. Appl. 27 (1994) pp. 59- 90. 

[9] M. Ida and T. Yabe, Implicit CIP (cubic interpolated 
propagation) method in one dimension, Comput. Phys. 
Commun. 92 (1995) pp. 21-26. 



[10] T. Utsumi, Differential algebraic hydrodynamics solver 
with cubic-polynomial interpolation, CFD J. 4 (1995) 
pp.225-238. 

[11] F. Xiao, T. Yabe and T. Ito, Constructing oscillation pre- 
venting scheme for advection equation by rational func- 
tion, Comput. Phys. Commun. 93 (1996) pp. 1-12. 

[12] F. Xiao, T. Yabe, G. Nizam and T. Ito, Constructing 
multi-dimensional oscillation preventing scheme for ad- 
vection equation by rational function, Comput. Phys. 
Commun. 94 (1996) pp.103-118. 

[13] M. Ida, An improved unified solver for compressible and 
incompressible fluids involving free surfaces. Part I. Con- 
vection, Comput. Phys. Commun. 132 (2000) pp.44-65. 

[14] M. Ida, "Non-smooth solution of density jump at inter- 
faces with discontinuous interpolation and a level set ap- 
proach", in: Proc. 10th Symp. on Comput. Fluid. Dy- 
nam., Tokyo, Japan, 1996, pp. 382-383 (in Japanese. Most 
all of this content is included in Ref . 113] ) . 

[15] M. Ida, "Free-surface flow simulation by an interpolation- 
extrapolation hybrid scheme", in: Proc. 10th Comput. 
Mech. Conf., Tokyo, Japan, 1997, pp. 17-18 (in Japanese. 
Most all of this content is included in Ref. fl3| ) . 

[16] T. Aoki, Comput. Phys. Commun. Interpolated differen- 
tial operator (IDO) scheme for solving partial differential 
equations, 102 (1997) pp. 132-146. 

[17] R. Tanaka, T. Nakamura and T. Yabe, Constructing 
exactly conservative scheme in anon-conservative form, 
Comput. Phys. Commun. 126 (2000) pp.232-243. 

[18] R. Tanaka, T. Nakamura and T. Yabe, "Constructing 
exactly conservative scheme in non-conservative form", 
in: Proc. 13th Symp. on Comput. Fluid. Dynam., Tokyo, 
Japan, 1999, pp. 1-5 (in Japanese, unpublished). 

[19] J. C. Clements, Convexity-preserving piecewise rational 
cubic interpolation, SIAM J. Numer. Anal. 27 (1990) 
pp.1016-1023. 



